Relaxation Time and Relaxation Function of Quark-Gluon Plasma with Lattice QCD 



OO 
O 

o 

(N 
>—> 



o 

=5 
S3 



OO 

o 

OO 

o 



X 



Kyosuke Tsumura 

J Analysis Technology Center, Fujifilm Corporation, Kanagawa 250-0193, Japan 

(Dated: June 9, 2008) 

We propose a method which enables a QCD-based calculation of a relaxation time for a dissipative 
current in the causal and dissipative hydrodynamic equation derived by Israel and Stewart. We 
point out that the Israel-Stewart equation is not unique as a causal and dissipative hydrodynamic 
equation, and the form of the causal and dissipative hydrodynamic equation is determined by the 
shape of a spectral function reflecting the properties of elementary excitations in the system we 
consider. 

Our method utilizes a relaxation function, which can be calculated from QCD using the linear 
response theory. We show that the relaxation function can be derived from a spectral function 
for a microscopic representation of the dissipative current. We also show that the Israel-Stewart 
equation is acceptable only as long as the calculated relaxation function is approximated well by 
a exponentially damping function, and the relaxation time can be obtained as its damping time 
constant. Taking a baryon-number dissipative current of a plasma consisting of charm quarks and 
gluons as a simple example, we present the first calculation of the relaxation function with use of the 
spectral function derived employing the quenched lattice QCD together with the maximum entropy 
method. The calculated relaxation function shows a strongly-oscillation damping behaviour due to 
the charmed vector hadron J/SP surviving above the deconfinement phase transition temperature in 
QCD. This result suggests that the applicability of the Israel-Stewart equation to the baryon-number 
dissipative current of the charm quark-gluon plasma is quite doubtful. We present an idea for the 
improvement of the Israel-Stewart equation by deriving the hydrodynamic equation consistent with 
the strongly-oscillation damping relaxation function. 

PACS numbers: 12.38.Mh, 47.75.+f, 12.38.Gc 



I. INTRODUCTION 

The time evolution of the hot matter of quarks and 
gluons created in the Relativistic Heavy Ion Collider 
(RHIC) experiments at Brookhaven National Laboratory 
is well described by the relativistic hydrodynamic equa- 
tion with almost no dissipative effects [H 0, Q . This fact 
means that the quark-gluon plasma (QGP) slightly above 
the deconfinement phase transition temperature Tq be- 
haves as an almost ideal fluid, of which viscosity and heat 
conductivity, i.e., transport coefficients, are negligiblely 
small. This interesting non-equilibrium property of QGP 
as an almost ideal fluid is considered to be one of the 
properties of a strongly interacting system, and consis- 
tent with hadronic excitations surviving in QGP shown 
by several groups using the lattice QCD 0, [Hi H 0] • 

In order to reach a full understanding of this non- 
equilibrium property of QGP, we need to study the RHIC 
experimental data using the relativistic hydrodynamic 
equation with finite, dissipative effects. 

At present, a good candidate for the relativistic hydro- 
dynamic equation with dissipative effects is considered to 
be the Israel-Steart equation [f| @ instead of the stan- 
dard equations by Eckart [l(| and Landau and Lifshitz 
plj . This is because the Israel-Stewart equation holds 
the causality in contrast with the Eckart and Landau- 
Lifshitz equations. 

To assure the causality, the Israel- Stewart equation 
contains higher derivative terms with relaxation times, 
i.e., new parameters other than transport coefficients, 



which are dependent on materials we consider. An impor- 
tant matter is to determine these relaxation times from 
microscopic dynamics governing the materials. 

In 1979, Israel and Stewart @ calculated the relax- 
ation times of a rarefied gas by using the relativistic 
Boltzmann equation with the Grad's fourteen moment 
approximation I12B . These quantities have been used by 
several groups [l3lll4llTH to carry out (2+l)-dimensional 
hydrodynamic simulations based on the Israel-Stewart 
equation. 

In 2008, Natsuume and Okamura [l(| derived the re- 
laxation times of a strongly interacting supersymmetric 
gauge-theory plasma by utilizing the dispersion relations 
of the hydrodynamic modes obtained with use of the 
AdS/CFT duality. However, a calculation of the relax- 
ation times of QGP based on QCD is unfortunately not 
yet done. 

The aim of this paper is (1) to construct a method 
which enables us to determine the relaxation times from 
the microscopic theory including QCD, and (2) to present 
the first lattice QCD calculation of the relaxation times 
of QGP on the basis of this method. 

To construct the method, we utilize a relaxation junc- 
tion, which gives a linear but non-local (non-Markovian) 
relation between an external field applied to the system 
and the current induced in the system [Tt], EH • The re- 
laxation functions contain full information about non- 
equilibrium dissipative processes, e.g., transport coeffi- 
cients and relaxation times. Moreover a microscopic cal- 
culation of the relaxation functions is straightforward us- 
ing the linear response theory. 
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In this paper, we first express the relaxation function 
in terms of a spectral function for a conserved current, 
which can be calculated from the microscopic dynamics 
by a quantum-field theoretic method. This expression 
tells us that the relaxation function is determined by a 
whole behaviour of the spectral function in contrast to 
the transport coefficient defined as the low-frequency and 
long-wavelength limit of the spectral function. We show 
that, only as long as the calculated relaxation function 
shows an exponentially damping behaviour, the damp- 
ing time constant can be identified as the relaxation time 
in the Israel-Stewart equation. Taking a baryon-numbcr 
dissipative current of a plasma consisting of charm quarks 
and gluons (referred to as the charm quark-gluon plasma 
later on) as a simple example, we calculate the relax- 
ation function from the spectral function obtained us- 
ing the lattice QCD together with the maximum entropy 
method (MEM) [19] . The obtained relaxation function 
shows a strongly-oscillation damping behaviour instead 
of an exponentially damping one, due to the charmed 
vector hadron J/ty surviving above Tq- This disagree- 
ment suggests that the applicability of the Israel-Stewart 
equation to the baryon-number dissipative current of the 
charm quark-gluon plasma is doubtful. We show an idea 
for the improvement of the Israel-Stewart equation by 
deriving the hydrodynamic equation consistent with the 
strongly-oscillation damping relaxation function. 

The conclusion of this paper is that the Israel-Stewart 
equation is not unique as the causal and dissipative hy- 
drodynamic equation, and the form of the causal and 
dissipative hydrodynamic equation applicable to QGP is 
determined by the shape of the spectral function reflect- 
ing the properties of hadronic excitations in QGP. This 
indicates an important link between the non-equilibrium 
properties of QGP and the hadron spectroscopy at finite 
temperature. 

This paper is organized as follows. In section [Til we 
construct a method which enables a microscopic calcu- 
lation of the relaxation times in the Israel-Stewart equa- 
tion. For this purpose, the relaxation functions are uti- 
lized. In section IIII1 we present the first calculation of 
the relaxation function for the baryon-number dissipa- 
tive current in QGP, with use of the Monte Carlo simu- 
lation based on the lattice QCD. We discuss whether the 
exponentially damping relaxation function of the Israel- 
Stewart equation can be justified from a microscopic 
point of view. Section HVl is devoted to the summary and 
the concluding remarks. In appendix [Al taking a shear 
viscous dissipative process as a typical example, we show 
the detailed derivation of the microscopic representation 
of the relaxation function. 



II. FORMULA FOR RELAXATION FUNCTION 

In this section, we give a method which enables us 
to determine the relaxation time from the microscopic 
theory including QCD. Our discussion is based on the 



relaxation function established in the phenomenological 
relaxation theory [13] • The microscopic representation 
of the relaxation function is given in the linear response 
theory [13, 111]. 



A. Non-Markovian Constituent Equation Derived 
from Israel-Stewart Equation 

To explain the physical meaning of the relaxation func- 
tions, let us start from the Israel-Stewart equation for a 
simple system, where the conserved quantities are the en- 
ergy and the momentum, and the dissipative effect comes 
solely from the shear viscosity. We denote the energy- 
momentum tensor by 



(1) 



where e = e(T) and p — p(T) are the internal en- 
ergy and the pressure at the temperature T, respec- 
tively, the flow velocity, 7r M1/ the shear viscous stress 
and = g^" — u v the projection operator with 

g^ v = diag(+l, — 1, — 1, — 1). Here w M is normalized as 
u^u^ = 1, and tt^ is satisfied with ir^ — tt^p and 
ir^ — = 0, which are reproduced by tt^ — 
A^ pa TT pa with the symmetric traceless projection op- 
erator A^P a = l/2(A^ ^ PA•" 7 + A^ Ia A l, P-2/3A^ n ^ , A P<T ) 
and the shear viscous stress before the projection 7r pcr . It 
is noted that the Israel-Stewart equation is the system 
of the nine differential equations with respect to the nine 
independent variables, T, u^ 1 and tt^. The equations 
consist of the continuity equation, 



d^ v = 0, 
and the relaxation equation [1, [9[ , 

7V DtT p<7 + TT pcr = 2 T] dpUa, 



(2) 



(3) 



where 77 and T n denote the shear viscosity and the re- 
laxation time of the shear viscous dissipative process, re- 
spectively, and D = = d / dr with the proper time 

T. 

The solution of Eq. ([3]) for the shear viscous stress reads 



TTpcr (t) = J 6.T 1 



r')/r 



2d p u a (T').(4) 



1 pa 



at the proper 



Equation (gj) tells us that vr^ = A^ 1 
time r depends on the history of the external field 2 d p u a 
from past to present, and the weighting function is given 
by ^-e-( r ~ r )/ T ". Notice that Eq.® is equivalent to 
that derived in the Stewart's non-local thermodynam- 
ics 1201. and has been recently rediscovered by Koide et 
al. [21(. The meaning of Eq.® becomes clearer from 
the viewpoint of the phenomenological relaxation the- 
ory [13] , where the linear but non-Markovian constituent 
equation, 



1 pa 



(r) 



J At 1 R{r-T , )2d p u r7 (T), 



(•5) 



3 



is well established and R(t) is identically the relaxation 
function. Therefore, we find that the Israel-Stewart 
causal hydrodynamics is equivalent to the special case 
of the phcnomcnological relaxation theory with 

R(r)=e(r)^-e- T ^, (6) 
i.e., the exponentially damping relaxation function. 



B. Microscopic Representation of Relaxation 
Function Based on Spectral Function 

We first notice that the linear response theory [It], Gil 
enables us to calculate the relaxation function directly 
from the microscopic dynamics. We now consider the 
relaxation function of a generic dissipative process, 



J(t, x) — Jdud 3 y R(t — u, x — y) X(u, y), 



(7) 



where X(t, x) denotes the external field, J(t, x) the cur- 
rent induced by X(t, x), and R(t, x) the relaxation func- 
tion. The linear response theory allows that the micro- 
scopic derivation of the relaxation function is straightfor- 
ward. In fact, as shown in appendix^ we can derive the 
following relation, 

R(t,x) = 6(t) [ ds8(t- s)i([J(t, x), J(s, 0)]), (8) 



where J(i, x) is a Heisenberg operator and gives the mi- 
croscopic representation for the induced macroscopic cur- 
rent J(t, x), and the bracket indicates the thermal aver- 
age. Notice that the integrand of the right hand side 
of Eq.® is the retarded Green function, which has the 
following spectral representation, 



9(t)i([J(t, x), J(0, 0)]) 
d4fc -ik°t+ik-* [ duJ Mf^ fc ) 



(27T) 



A: 



IE 



(9) 



where A(u>, fc) is the spectral function and e the infinites- 
imal positive constant. Equations ([8]) and ([9]) lead us to 
the relation between R(t, x) and A(u>, fc) given by 



R(t, x) = 2 9(t) 



dwd 3 fc 



Here we notice the two important properties concerning 
the complex admittance R(k°, fc): (i) The low- frequency 
and long-wavelength limit of the complex admittance, 
R(0, 0) = itA(l:, 0)/u/| , gives the transport coeffi- 
cients. In fact, with use of this representation, many 
groups [H, HI, 0, HE] have calculated the transport co- 
efficients of QGP employing the lattice QCD. (ii) Using 
i?(fc°, fc), we can write down the generic relaxation equa- 
tion equivalent to the non-Markovian constituent equa- 
tion ©, 



R-^k = +id t , fc = -iV)J(t, x) = X(t, x). 



(12) 



By combining Eqs. pTj) and (H2), we can construct the 
relaxation equation from a given spectral function. 

The above observation tells us the following three 
points: (a) The relaxation function in QGP can be cal- 
culated from the spectral function for the microscopic 
current of QCD. (b) The Israel-Stewart equation can be 
valid for the description of QGP only as long as the relax- 
ation function calculated in the microscopic way is an ex- 
ponentially damping function, and the relaxation time of 
QGP can be obtained as its damping time constant, (c) If 
the relaxation function calculated from QCD is different 
from that derived from the Israel-Stewart equation, we 
should carry out the causal hydrodynamic calculation by 
incorporating the relaxation equation (|12p or the non- 
Markovian constituent equation (JT]) into the continuity 
equation. 



C. Spectral Function Consistent with 
Israel-Stewart Equation 

It is interesting to investigate which spectral functions 
reproduce the exponentially damping relaxation function 
in the Israel-Stewart equation. Taking the shear viscous 
dissipative process at the rest frame as a typical example, 
we find it to be a Lorentzian function, 



-A(w, fc) = - — V^- 



(13) 



In fact, it can be checked that Eqs. (fT0|) and (( 13|) give the 
exponentially damping relaxation function, 



e 



— iut-\-ik-x 



■ A(w, fc). (10) & n d the correspondent complex admittance reads 



This formula assures that the relaxation function can be 
obtained as the Fourier transformation of ttA(oj, k)/u>, 
i.e., the real part of the complex admittance [17| . which 
is defined by 



R(k°, fc) = 



n 



l-ik°- 



(14) 



(15) 



It is noteworthy that Eq. (|12p with this complex admit- 
tance reproduces the Israel-Stewart-type relaxation equa- 
tion, 



R{k°,k) = Jdtd 3 xR(t, x)e lk0t 

-i A(uj, fc) 



(l + 7v<9 t )J(i, x)=r)X(t, x). 



(16) 



duj ■ 



— k° — is 



(11) 



Here we replace dt with D = u^d^ — d/dr in order to 
obtain the equation at an arbitrary frame from that at 
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the rest frame. Therefore it is found that the ansatz of 
the Lorentzian-type spectral function shown in Eq. (|13p 
is an essential point of the Israel-Stewart equation. No- 
tice that there is no reason to believe that the spectral 
function for the current of QCD is approximated well by 
the Lorentzian function in the low-frequency and long- 
wavelength region. 



III. RESULTS FROM LATTICE QCD 

In this section, we demonstrate the first lattice QCD 
calculation of the relaxation function. From the view- 
point of availability of lattice data, we choose the relax- 
ation function for the baryon-number dissipative process. 
The continuity equation reads 

d„N» = d„(nu> t + iS') = 0, (17) 
and the non-Markovian constituent equation 



v»(t, x) 



A^(i, x) Jdud 3 yR{t 



y 



xdApB/T)(u, y), (18) 



where n and /j,b denote the baryon-number density and 
the baryon-number chemical potential, respectively. This 
relaxation function can be calculated from the spectral 
function of the spatial components of the baryon-number 
current, = J2fZ\ X) c =i V'/cT^V'/c where ip fc de- 
notes the /-flavour c-colour quark field, 7 p (/i = 0, 1, 2, 3) 
the gamma matrix, 'ipfc = V'/cT : -^f the number of 
the flavour and Nc(— 3) the number of the colour, re- 
spectively. As a simple example, we consider the single- 
flavour deconfined system of which constituents are glu- 
ons and charm quarks (N-p — 1), i.e., the charm quark- 
gluon plasma. Furthermore we treat a spatially inte- 
grated relaxation function, 



R(t) = Jd 3 xR(t, x) = 2 0(t) J- 



'fe-^AM, (19) 

Z7T UJ 



where A(u>) = A(u>, 0) is a zero-momentum projected 
spectral function. This relation indicates that the calcu- 
lation of R(t) requires an explicit form of A(u). Therefore 
we need to reconstruct A(u>) from the temporal correla- 
tion function, 



3 

D(t) = 1/3 J2 J d3x ( Jl (-^ x ) J X0, 0) ), 



(20) 



obtained using Monte Carlo simulation based on the lat- 
tice QCD. To get the full information about A(u>), we 
utilize MEM ml. 



A. Outline of MEM Analysis 

Here we shall give a brief account of MEM. The tem- 
poral correlation function D(t) and the associated A(ui) 



are related by the Laplace transformation, 

/>oo 

D(t)= &uK{t,uj)A{u)), (21) 

JO 

where K(t, to) = cosh[w (r-l/2T)]/sinh(cj/2T) is a ker- 
nel function. Since the number of data points in the 
temporal direction is finite, the direct inverse Laplace 
transformation from D(t) to A{u>) is an ill-posed prob- 
lem. MEM is a method to overcome this difficulty on 
the basis of the Bayes' theorem in the theory of statis- 
tical inference. Using MEM, we can infer A(ui) as the 
most probable one which achieves a balance between re- 
producing the lattice data D{t) and being kept not so 
far from the default model to(w), defined by a plausible 
form of A(lu). MEM has proved to be a powerful tool to 
inve stig ate A(ui) of mesons and baryons at both of T = 

[HiEg and T + fliJLIS- 

Recently, the author [27| has improved MEM to avoid 
a numerical instability present in the study of the shear 
viscosity of the gluon plasma. An essential point of the 
improved MEM analysis is that the Laplace transforma- 
tion combining D(t) and A(uj) is converted as follows, 



D(r) 



dw 



K(t, u) - 

IT 



-A(u) 



(22) 



Regarding u> K(t, uj)/tt as a new kernel function, we can 
infer the real part of the complex admittance 7r A(u>)/u> 
directly from D(r). Further, in the study of the scalar 
tetraquark 28], the extension of this method, such as 
D{t) = f^duj [K{t, w)A(w)] [A(w)/A(w)] with A(w) = 
uj s , has also been used to reduce the uncertainty of the 
inferred A(u>). Aarts et al. [25| independently have pre- 
pared this method, and have applied this method to cal- 
culate the electrical conductivity of QGP. 

In FigfTJ we show the schematic overview which sum- 
marizes the way to obtain the relaxation function and the 
relaxation time of QGP on the basis of the lattice QCD 
and the improved MEM analysis. 

Then we explain the default model used in the present 
analysis. The default model m(u>) is usually given as 
the form predicted by the perturbativc QCD calculation 
valid in the high-energy region, lu — > oo. Due to the tree- 
level result, A(lo) — moui 2 with mo = l/4n 2 , the default 
model reads 



m{ui) 



(23) 



Such a monotonically increasing default model seems 
suitable to avoid the artificial peaks due to the choice of 
the default model, and using m mono (aj) we can identify 
any peak structure of the inferred A(u>) as that derived 
from the lattice data not the default model. In this paper, 
we are interested in the applicability of the Israel-Stewart 
equation to QGP, i.e., the validity of the ansatz that the 
spectral function for the conserved currents is approxi- 
mated by a Lorentzian function shown in Eq. (|13[) . To 
test this ansatz, we carry out the MEM analysis based 
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Relaxation Function R(t) ] I S Relaxation Time X 

in Israel-Stewart equation 



t 



Formula in Eq.(lO) 



Spectral Function irA(a))/(B 




t 



MEM based on Eq.(22) 



if R(t) oc e - t/T 
for t > 



Correlation Function D(t) 



FIG. 1: The schematic overview to calculate the relaxation 
times in the Israel-Stewart causal and dissipative hydrody- 
namic equation employing the lattice QCD together with the 
maximum entropy method (MEM). The correlation function 
can be calculated in lattice. 



on not only the monotonically increasing default model 
m mono (w) but also the following Lorentzian-type one es- 
pecially in the low-energy region, 



m(to) = m\ 



M - \ ( 



= — ( 1 — tanh 



LU — LUq 



X 



1 + tanh 



Auo ) 1 



A 



it mo w,(24) 



where lo and Aw denote the center value and width of 
the border between the strongly interacting low-energy 
region and the asymptotically free high-energy one, re- 
spectively, and x an d T J denote the transport coefficient 
and relaxation time of the baryon-number dissipative 
process, respectively. If it A{lo)/lo inferred by the MEM 
analysis based on m(u>) = mi OIe (cj) is quite different from 
the Lorentzian function in the low-frequency region, the 
applicability of the Israel-Stewart equation to QGP is 
doubtful. 



B. Numerical Result and Discussion 
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FIG. 2: (a) The real part of the complex admittance 
7T A(lo)/loT inferred using the improved MEM analysis based 
on m{uo) = m mono (cj) in Eq, (l23fl . (b) The result of the im- 
proved MEM analysis based on m(uo) — mi ore (u;) in Eq. (|24p 
with loo = 1.0 GeV, Alo = 2.0 GeV, x = 5.0 GeV and 
r.j = 0.1 GeV -1 . The dash line denotes the default model 
m(u))/T and the solid line the inferred it A(lo)/loT. The hori- 
zontal position and length of the bars indicate the region over 
which ttA(lo)/loT is averaged, while the vertical height of the 
bars denotes the standard deviation in the averaged value of 
ttA(lo)/loT in the interval. 



The system considered here is the charm quark-ghton 
plasma at T > Tc- We use the quenched approximation 
in calculating D(t). A detailed set up of the lattice QCD 
numerical simulation is identical to that in Ref.[3]. We 
use the naive plaquette gauge action and the standard 
Wilson quark action with the gauge coupling constant 
f3 = 7.0, the bare anisotropy (q = 3.5 and the fermion 
anisotropy 7f = k t /k ct = 3.476 with the spatial (tempo- 
ral) hopping parameter n a = 0.08285 (k t ). These param- 
eters reproduce the anisotropic lattice a a /a T = 4 with the 
temporal (spatial) lattice spacing a T = 9.75 x 10~ 3 fm 
(a a ) in the quenched-level simulation Q. The adopted 
lattice size is 20 3 x 46, where the spatial size L a — 0.78 fm 
and the temperature T = 1.62 Tc- This spatial size seems 
small. By varying the spatial boundary condition, how- 
ever, we have checked that L a is large enough to simulate 



the charm quark-gluon plasma in lattice Q ■ The number 
of the gauge configurations is 100. 

The improved MEM analysis is done using the lat- 
tice data D(t) (r = Ti,--- , r/) with Ti/a T — 4 and 
Tf/a-r — 42 and the default models m mono (uj) in Eq. (|23|) 
and mi orc (a;) in Eq.flMJ). In Fig[21 we show the inferred 
7r A(u>) /o/Sj which have been checked to be stable against 
the change of (Tj,r/) as (5,41), (6,40) and (7,39). Figure 
[2]-(a) indicates that it A(lo)/uj has no strength in the low- 
frequency region lo < 3 GeV and the peak structure with 
large width around u ~ 4.5 GeV. Notice that due to 
m mono (u)) having no strength in lo < 3GeV and taking a 
zero value at to = 0, we cannot get any information about 
the applicability of the Israel-Stewart equation to the 
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baryon-number dissipative process of the charm quark- 
gluon plasma and the smallness of the transport coef- 
ficient of the baryon-number dissipative process x from 
this result. On the other hand, the peak structure with 
large width around u> ~ 4.5 GeV can be considered to 
be a signal of the charmed vector hadron J/\I/ surviving 
above Tc- This result is in good agreement with those 
of the several groups 0, H, H, 0| ■ Figure [21- (b) shows that 
7r A(u)/cu is quite different from mi mc (u), the Lorentzian 
function in the low-frequency region. This result sug- 
gests that the applicability of the Israel-Stewart equation 
to the baryon-number dissipative process of the charm 
quark-gluon plasma is questionable. We have checked 
the stability of the inferred it A{uS)/u> against the change 
of wo, Aw, x an d Tj. 

We calculate the relaxation function expected to be 
quite different from the exponentially damping function. 
For this purpose, we carry out the Fourier transformation 
of the inferred it A(lo)/uj in Fig[2]-(a). Since it A(lo)/lo ex- 
hibits a linear divergence of to for to — > oo, the numerical 
Fourier transformation is unstable. To avoid this numer- 
ical instability, we convert Eq. (fT9| into 




-100 



0.00 0.05 0.10 0.15 0.20 0.25 0.30 

t[fm] 



FIG. 3: The solid line denotes the relaxation function 
R(t)/T 2 with T = 1.62 T c . The meaning of the horizontal 
position and length of the bars is explained in the caption 
of Fig(2] R(t) /T 2 exhibits a strongly-oscillation damping be- 
haviour within the error bars, instead of the exponentially 
damping behaviour derived from the Israel-Stewart equation. 



R(t) = 2 6{t) 



mo 
t 2 



dui cos tot 



-A(u) 

10 



(25) 



where we have used A(to) 
relation, 



dwe" Mt ( 



-A{— to) and the following 



for t f 0. 



(26) 



In Eq. ([2"5]l . Lo cut denotes the cut-off frequency de- 
fined as ttA(lo)/lo — m mon o(w) ~ for to > Lo cut . 
Here we use ui cu t = 100 GeV > n/a T — 63 GeV. 
In Fig[3l we show the calculated R(t) together with 
the error bars. The singularity of R(t — 0) is in- 
evitable because it is originated from the divergence of 
ttA(lo)/lo at to — > oo due to the composite operator 

J»{-iT,x) = T,fZiT,?=i$fc(-iT,x)j»ip(-iT,x) fc 
being ill-defined at the same coordinate. Therefore, we 
concentrate on R(t) for t > a T — 9.75 x 10~ 3 fm. As 
expected, R(t) is not a exponential damping function 
but a strongly-oscillation damping function. This oscilla- 
tion damping behaviour of R(t) is attributed to the peak 
structure with the finite width in ttA(lo)/lo correspon- 
dent to the surviving J/^S at T — 1.62 Tc- The relaxation 
time Tosc of the oscillation damping function is approxi- 
mated by the inverse of the width of the low-lying peak 
structure. Accordingly, we can estimate r osc ~ 0.1 fm 
from Fig[21-(a), which seems to be consistent with R(t) 
in FIG. [31 The disappearance of J/vf at the higher tem- 
perature would drive the drastic change of R(t) and t osc . 



C. Relaxation Equation Applicable to Strongly 
Interacting System 

Finally we discuss the form of a relaxation equation 
corresponding to the spectral function with the peak 
structures in the finite-frequency region. Here we con- 
sider the spectral function for the shear viscous dissipa- 
tive process at the rest frame as a typical example, 



■A(u, k) 



r 7 (^+ r - 2 )/2 , r?(^+r- 2 )/2 



(iO + LO^) 2 + TV 



=2>(27) 



where w T (> 0) denotes the peak position. Notice that 
Eq.([27| agrees with Eq.(l3]) by the setting of uj^ = 0. 
Using Eqs. lfTTl) and ([27]) . we have the correspondent com- 
plex admittance, 



CO 2 - + T- 2 



-ik° 



j 2 + (-ik° + t^ 1 ) 



(28) 



By combining Eqs.([T2j) and ([28|) . we find the form of the 
relaxation equation corresponding to the spectral func- 
tion with the two peak structures at to = iav to be 



2 2 
,2 _2 



+ (1 + TV d t f J(t, X) 



rj(ivi4 + l)(l + T^d t )X(t, x). 



(29) 



It is found that the temporal derivative dt appears also 
in the right hand side of Eq. (|2"9")) , in contrast with the 
Israel-Stewart- type relaxation equation ^rJ]) . We identify 
Eq. ([29l) as a more probable form of the relaxation equa- 
tion applicable to the strongly interacting system includ- 
ing QGP slightly above Tc, of which spectral function 
has peak structures in the finite-frequency region. 
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As shown in subsection IIIIB1 the baryon- number dis- 
sipative process of the charm quark-gluon plasma cannot 
be described by the Israel-Steart equation. We should 
carry out the causal hydrodynamic simulations by incor- 
porating the relaxation equation based on Eq. (|29p not 
Eq. (fT6|) into the continuity equation (fT7|) . 



IV. SUMMARY AND CONCLUDING 
REMARKS 

First, in this paper, we have pointed out that the 
Israel-Stewart causal hydrodynamic equation is accept- 
able only as long as the relaxation function is a exponen- 
tially damping function, and the relaxation time can be 
obtained as its damping time constant. Using the linear 
response theory, we have derived the formula where the 
relaxation function is represented by the spectral func- 
tion for the conserved current, which can be calculated 
from the microscopic dynamics. Furthermore, we have 
discussed the systematic way to construct a causal hy- 
drodynamic equation corresponding to the calculated re- 
laxation function. 

Then, taking the baryon-number dissipative process of 
the charm quark-gluon plasma for a simple example, we 
have calculated the relaxation function from the spectral 
function obtained using the quenched lattice QCD to- 
gether with the maximum entropy method. The obtained 
relaxation function shows the strongly-oscillation damp- 
ing behaviour within the error bars, not the exponential 
damping. The origin of the strongly-oscillation damping 
behaviour is the charmed vector hadron J/^S surviving 
at T — 1.62 Tc- This fact suggests that the applicabil- 
ity of the Israel-Stewart equation to the baryon-number 
dissipative process of the charm quark-gluon plasma is 
questionable. We have shown an idea for the improve- 
ment of the Israel-Stewart equation by deriving the re- 
laxation equation consistent with the strongly-oscillation 
damping relaxation function. 

We conclude that the Israel-Stewart equation is not 
unique as the causal and dissipative hydrodynamic equa- 
tion. By calculating the relaxation function from the 
spectral function with use of the formula in Eq. (TII)]) , 
we can investigate the applicability of the Israel-Stewart 
equation. It is noteworthy that the form of the causal 
hydrodynamic equation with dissipative effects applica- 
ble to QGP is determined by the shape of the spectral 
function reflecting the hadronic excitations in QGP. This 
indicates an important link between the non-equilibrium 
properties of QGP and the hadron spectroscopy at the 
finite temperature. 

In the next paper [30[, we will investigate the prop- 
erties of the relaxation equation corresponding to the 
spectral function with the peak structures in the finite- 
frequency region. It seems an interesting subject to de- 
rive the relativistic hydrodynamic equation by combining 
the relaxation equation shown in Eq. (|29p and the con- 
tinuity equation based on the energy-flow equation by 



Landa u [111 or the particle-flow equation by Tsumura et 
al. [3lL |32| ]. which have a stable equilibrium solution. 

To discuss the more realistic system than the charm 
quark-gluon plasma, we must investigate the contribu- 
tion of up, down and strange quarks to the spectral 
function for the baryon-number conserved current. The 
light vector hadrons composed of these quarks slightly 
above Tc seem to play a more essential role for the slow 
and long- wavelength fluid behaviour of QGP than J/ty. 
The answer to this conjecture will be reported elsewhere. 
In addition, the calculation of the relaxation function 
and/or the spectral function for the shear and bulk vis- 
cous dissipative process of QGP is an interesting future 
study using the lattice QCD with the maximum entropy 
method or another non-perturbative method, such as the 
QCD sum rule or the AdS/CFT duality. 
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APPENDIX A: MICROSCOPIC 
REPRESENTATION OF RELAXATION 
FUNCTION 

In this section, we show the derivation of Eq. ((SJ) based 
on the linear response theory [ijj together with the 
non-equilibrium statistical operator method proposed by 
Zubarev 

We consider a statistical operator p(t) = | ^(t) ){ i S{t) |, 
where | ^(t) } denotes a state of the system. The time 
evolution of p(t) is governed by the Liouville equation, 

d t p{i)+i[H,p{i)\=V, (Al) 

or the Neumann equation, i.e., the equation converted 
from Eq. fATj) . 

d t \np(t)+i[H,\np(t)}=0, (A2) 

both of which are equivalent to the Schrodinger equation, 
id t \V(t)) = H \%{t)) with H being the Hamiltonian. It 
is noted that Eqs. (|Al[) and (|A2|) hold the time-reversal 
symmetry. Therefore, the entropy of the system defined 
by S(t) = Tr p(t) In p(t) is normally conserved. 

To obtain the non-equilibrium statistical operator 
Pneq{t), which exhibits the increase of the entropy, the 
breaking of the time-reversal symmetry is needed. In 
Ref. [l|| , Zubarev has proposed to construct p neq (t) as the 
solution of the Neumann equation (|A2|) under the time- 
reversal asymmetric boundary condition p{t — > —00) = 
Pqeq{t — > — °o|0), where p qoq (t\0) denotes the quasi- 
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equilibrium statistical operator defined by 



Pqeq(tl\t 2 ) OC exp 



^xa^h, x)f"°(t 2 , x) 



-(A3) 



with Trp qeq (ii|t 2 ) = 1. Here, T pu (t, x) denotes a 
energy-momentum tensor operator in the Heisenberg pic- 
ture, which agrees with that of the Schrodinder picture 
at t = 0. a^t, x) denotes a Lagrange multiplier charac- 
terizing the system. The physical meaning of a M (i, x) 
becomes clear from the setting of a Al (i, x) — g^o/T, 
which reproduces the equilibrium statistical operator 
p cq ee e~ d / T /Tre-^/ T with H = Jd 3 xf 00 (0, x). 

It is noted that imposing the time-reversal asymmet- 
ric boundary condition p(t — »■ — oo) = Pq C q (t -oo|0) 
to Eq. (|A2p is equivalent to adding a time- reversal asym- 
metric source term to the right hand side of Eq. (|A2j) , 



8 t hxp{t) +i[H, \np(t)} = -e In \p(t) / p qcq (t\0)\ , (A4) 

where e denotes an infinitesimal positive constant and 
the limit of e — > +0 follows the thcrmodynamical limit. 
Pneq{t) can be obtained as /5 ncq (i) = lim e ^ + o Pe(t), where 
p e (t) denotes the solution of Eq. (|A4j) and lim e ^ + o is sup- 
pressed later on. The solution of Eq. (|A4[) reads 



lnpneq(i) 



dse 



e(a-t) 



ln/3 qeq (s|s - t). (A5) 



A partial integration of the right hand side of Eq. (|A5p 
leads us to 



Pn eq (t)=c A+B /Tre^ 



(A6) 



where 



A=- jd 3 xa ll (t, x)f»°(0, x), (A7) 
B = f dse< s -^ fcfxd^is, x)f pv {s-t 1 x). (A8) 



Here we have used <9 M T A " / = 0. 

Notice that B contains c^ t a„, i.e., the gradient of the 
Lagrange multiplier characterizing the system. We shall 
consider the system not so far from the equilibrium state, 
where d^a^, is small. It is important that under this 
consideration we can regard B as a perturbation term 
to A and utilize the linear response theory, where d^a v 
is identified as an external field. We expand p neq (t) up 
to the first order of B as follows, 



Pneq{t) = 1 



dA 



(B-(B) A )e- XA )p A ,(A9) 



where pa = e A / Tr e A and ( O } A = Tr pA O with an 
arbitrary operator O. Notice that p A = Pqcq(i|0). 



Using Eq. (|A9[) , we calculate the energy- momentum 
tensor realized in the non-equilibrium state, 



T^(t, x)=Tr Pneq (t)f^(0, x). 



(A10) 



We decompose T^ v (t, x) into the unperturbative term 
and the perturbative one as T pjV (t 1 x) = T^ v (t, x) + 
ST^(t, x), where 

T^(t,x) = (f^(0,x)) A , (All) 

ST^(t,x)= Jdud 3 yR pvpa {t, x\u, y) 

x d p a a {u, y). (A12) 

Here T^ v represents the perfect-fluid part of the energy- 
momentum tensor, while ST^ the dissipative part. 
R^ vpa {t, x\u, y) in Eq. (|A12[ ) is a relaxation function de- 
fined by 

R pupa {t, x\u, y) = 6{t - u) e - e (*-") 



dA(T^(0, x)e +XA T pa {u~t, y)e- AA ) A 



-XA > 



-(T»"(0, x)) A (T p °(u-t, y)) A 



(A13) 



Equation (j A13|) is one of the main results in this section. 

As the thermal average used in Eq. (IA13p is defined 
with use oi p A — /5qcq(^|0) not p eq , Eq. (|A13| ) seems incon- 
venient to calculate R pvp<J (t, x\u, y) employing a tech- 
nique established in the finite-temperature field theory. 
To derive R pvp<T (t, x\u, y) written by p cq , we shall make 
an approximation for Eq. (|A12| l. We consider a situa- 
tion where a fluctuation of a^t, x) is small. Notice that 
this approximation differs from that used in the linear 
response theory, and is implemented in a systematic way 
using the expansion, 



a M (i, x) — a^+e 8a M (t, x), 



(A14) 



where e is an infinitesimal constant different from e. We 
evaluate a leading order of ST^^, x) with respect to e. 
Due to dpa a (u, y) being the first order of e, the zeroth- 
order evaluation of R Pll,f " J (t, x\u, y) is sufficient. 
First, using a^it, x) = a^, we have 



(A15) 



where the energy-momentum vector operator = 
Jd 3 x r M °(0, x). With use of the unitary transformation 
for the temporal and spatial translation by P p , 

c +XA 6{t, x) e- XA = 6{t + iXa° 1 x + iXa), (A16) 

and the cluster property of the correlation function, 

( 6(0, x) 6(t, y)) A = { 6(0, x) ) A ( 6(t, y) ) A 

for |i| -> oo, (A17) 
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we can convert the right hand side of Eq. (|A13[) into 
0(t - u) e - £{t - u) f dX f ds(f" u (0,x) 



x d s f pa {u -t+(iX + s)a°,y + (iX + s) a) ) A .(A18) 

By replacing d s — ► —id\, integrating with respect to A 
and using the translational invariance, we find Eq. (|A18p 
to be 



9(t - u) e -<*-^ / ds 

J OO 

x i { [f p "(t -u,x-y), f pa {s a°, s a)] ) A . (A19) 



Using Eq. (|A19|) , we can rewrite Eq. (|A12[) as a Lorentz- 
covariant form, 

ST^(x) = Jd 4 yR p ^(x-y)d p a a (y), (A20) 
R pvpa {x) = 0(x°)e- ex ° / ds 

J — GO 

xi([f^(x),f^(sa)]) A . (A21) 

Then, using a p = u p / T (u^u p — 1), we derive 
R pupa (x) written by p eq . For this purpose, we utilize 
a decomposition of 5T pi/ (x) based on the irreducible ten- 
sors u p , A pv = g pv - u p u v and A pvpa = 1/2 ( A pp A va + 
A pa A vp ~ 2/3 A pu A pa ), namely, 

5T pv {x) = E{x) u p u v - U{x) A pv 
+Q p {x)u v + Q»{x)u p +ir p »(x), (A22) 

where E(x) = 8T ab {x)u a u b , IL(x) = -ST ab (x) A ab / 3, 
Q^x) = ST ab (x)u a A bf , and ir^(x) = 8T ab {x)A ab ^. 



Here we focus the shear viscous stress tt i11/ (x). Using 
Curie's principle (29^ and the spatial isotropy, we have 

{[^{x),fr°{su/T)]) A 
= A pupa -L ( [n« b (x) , n ab (su/T)} ) A , (A23) 

where tt^ = A^ ab f ab . This relation leads us to 

n^(x) = A M , pCT Jd 4 yR(x-y)2d p u°(y), (A24) 



0\ „-ex" 



R{x) = 6»(a; u )e 



1 



ds 



xi— ([7^(1), V(*«/r)]> A (A25) 



It is noted that, due to the integrand of R(x) in Eq. (|A25| ) 
being a Lorentz-scalar function, the integrand takes the 
same value for an arbitrary u p . Therefore, we can set 
= (1, 0, 0, 0) in Eq. (|A25|) without loss of generality, 
and using the spatial isotropy we arrive at 

R(x) = R(t, x) 



9(t)e 



x 6{t - s)i ( [T 12 (t, x) , Ti 2 (s, 0)] },(A26) 
where (O) = Trp cq O. Notice that the integrand of 
Eq. (|A26p is identical to the retarded Green function 
which has the spectral representation as shown in Eq. ([9]) . 
Equation (|A26|) is the other of the main results in this 
section. 
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